Introduction à la modélisation statistique bayésienne

Cours n°01 : Introduction à l'inférence bayésienne
Cours n°02 : Modèle Beta-Binomial
Cours n°03 : Introduction à brms, modèle de régression linéaire
Cours n°04 : Modèle de régression linéaire (suite)
Cours n°05 : Markov Chain Monte Carlo
Cours n°06 : Modèle linéaire généralisé
Cours n°07 : Comparaison de modèles
Cours n°08 : Modèles multi-niveaux
Cours n°09 : Modèles multi-niveaux généralisés
Cours n°10 : Data Hackathon

Le package brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le package lme4.

Reaction ~ Days + (1 + Days | Subject)

La partie gauche représente notre variable dépendante (ou outcome, i.e., ce qu'on essaye de prédire). Le package brms permet également de fitter des modèles multivariés (plusieurs outcomes) en les combinant avec c():

c(Reaction, Memory) ~ Days + (1 + Days | Subject)

La partie droite permet de définir les prédicteurs. L'intercept est généralement implicite, de sorte que les deux écritures ci-dessous sont équivalentes.

c(Reaction, Memory) ~ Days + (1 + Days | Subject)
c(Reaction, Memory) ~ 1 + Days + (1 + Days | Subject)

Si l'on veut fitter un modèle sans intercept (why not), il faut le spécifier explicitement comme ci-dessous.

c(Reaction, Memory) ~ 0 + Days + (1 + Days | Subject)

La première partie de la partie droite de la formule représente les effets constants (effets fixes), tandis que la seconde partie (entre parenthèses) représente les effets variables (effets aléatoires).

c(Reaction, Memory) ~ Days + (1 | Subject)
c(Reaction, Memory) ~ Days + (Days | Subject)

Le premier modèle ci-dessus contient seulement un intercept variable, qui varie par Subject. Le deuxième modèle contient également un intercept variable, mais aussi une pente variable pour l'effet de Days.

Lorsqu'on inclut plusieurs effets variables (e.g., un intercept et une pente variables), brms postule qu'on souhaite aussi estimer la corrélation entre ces deux effets. Dans le cas contraire, on peut supprimer cette corrélation (i.e., la fixer à 0) en utilisant ||.

c(Reaction, Memory) ~ Days + (1 + Days || Subject)

Les modèles précédents postulaient une fonction de vraisemblance Gaussienne. Ce postulat peut être changé facilement en spécifiant la fonction de vraisemblance souhaitée via l'argument family.

brm(Reaction ~ 1 + Days + (1 + Days | Subject), family = lognormal() )

Travailler avec des sujets humains implique un minimum de coopération réciproque. Mais ce n'est pas toujours le cas. Une partie non-négligeable des étudiants qui s'inscrivent pour passer des expériences de psychologie ne se présentent pas le jour prévu… On a voulu estimer la probabilité de présence d'un étudiant inscrit en fonction de l'envoi (ou non) d'un mail de rappel (cet exemple est présenté en détails dans deux blogposts, accessibles ici, et ici).

library(tidyverse)

data <- read.csv("data/absenteeism.csv")
data %>% sample_frac %>% head(10)

   reminder researcher presence absence total
1         1          4       92       3    95
2        -1         10       23      58    81
3        -1          4       61      34    95
4         1          6       82       6    88
5         1         10       76       5    81
6        -1          8       38      54    92
7         1          7       64      16    80
8        -1          1       16      86   102
9         1          5       65      18    83
10       -1          9       31      65    96

\[ \begin{aligned} y_{i} &\sim \mathrm{Binomial}(n_{i}, p_{i}) \\ \text{logit}(p_{i}) &= \alpha_{\text{researcher}_{[i]}} + \beta_{\text{researcher}_{[i]}} \times \text{reminder}_{i} \\ \begin{bmatrix} \alpha_{\text{researcher}} \\ \beta_{\text{researcher}} \\ \end{bmatrix} &\sim \mathrm{MVNormal}\bigg(\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \textbf{S}\bigg) \\ \textbf{S} &= \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \\ \end{pmatrix} \textbf{R} \begin{pmatrix} \sigma_{\alpha} & 0 \\ 0 & \sigma_{\beta} \\ \end{pmatrix} \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \beta &\sim \mathrm{Normal}(0, 1) \\ (\sigma_{\alpha}, \sigma_{\beta}) &\sim \mathrm{HalfCauchy}(0, 1) \\ \textbf{R} &\sim \mathrm{LKJcorr}(2) \\ \end{aligned} \]

Il s'agit du même modèle de régression logistique vu au Cours n°06, avec une fonction de lien logit, mais cette fois-ci sur plusieurs niveaux.

prior1 <- c(
    prior(normal(0, 1), class = Intercept, coef = ""),
    prior(normal(0, 1), class = b),
    prior(cauchy(0, 1), class = sd),
    prior(lkj(2), class = cor)
    )

mod1 <- brm(
    presence | trials(total) ~ 1 + reminder + (1 + reminder | researcher), 
    family = binomial(link = "logit"),
    prior = prior1,
    data = data,
    sample_prior = TRUE,
    warmup = 2000, iter = 1e4,
    chains = 4, cores = parallel::detectCores(),
    control = list(adapt_delta = 0.95),
    backend = "cmdstanr"
    )

mod1 %>%
    plot(
        combo = c("dens_overlay", "trace"), pars = c("^b_", "^cor_"), widths = c(1, 1.5),
        theme = theme_bw(base_size = 16, base_family = "Open Sans")
        )

Attention, les estimations ci-dessous sont dans l'espace log-odds…

posterior_summary(mod1, pars = c("^b_", "^sd_") )

                          Estimate Est.Error      Q2.5     Q97.5
b_Intercept              0.7892803 0.2632695 0.2518506 1.3060062
b_reminder               1.5623143 0.1649351 1.2244800 1.8816810
sd_researcher__Intercept 0.7950345 0.2320246 0.4642587 1.3598370
sd_researcher__reminder  0.4509342 0.1568272 0.2203203 0.8237991

Afin de pouvoir les interpréter il faut appliquer la transformation logit-inverse. Par exemple, la probabilité de présence en moyenne (i.e., quel que soit le chercheur et pour toutes conditions confondues) est égale à \( p = \exp(\alpha) / (1 + \exp(\alpha) ) \).

a <- fixef(mod1)[1] # extract the intercept
exp(a) / (1 + exp(a) ) # equivalent to plogis(a)

[1] 0.6876768

On s'est ensuite interrogé sur l'effet du mail de rappel. Ici encore, on ne peut pas interpréter la pente directement… mais on sait que \( \text{exp}(\beta) \) nous donne un odds ratio (i.e., un rapport de cotes).

fixef(mod1)[2, c(1, 3, 4)] %>% exp

Estimate     Q2.5    Q97.5 
4.769847 3.402396 6.564531 

Ce ratio est environ égal à 4.77. On peut l'interpréter en disant qu'envoyer un mail de rappel multiplie par 4.77 la cote (i.e., le rapport \( \frac{\Pr(\text{présent})}{\Pr(\text{absent})} \)).

Une manière de représenter les prédictions du modèle est de plotter directement quelques échantillons issus de la distribution a posteriori. On appelle ce genre de plot un spaghetti plot.

data %>%
    group_by(researcher, total) %>%
    data_grid(reminder = seq_range(reminder, n = 1e2) ) %>%
    add_fitted_samples(mod1, newdata = ., n = 100, scale = "linear") %>%
    mutate(estimate = plogis(estimate) ) %>%
    ggplot(aes(x = reminder, y = estimate, group = .iteration) ) +
    geom_hline(yintercept = 0.5, lty = 2) +
    geom_line(aes(y = estimate, group = .iteration), size = 0.5, alpha = 0.1) +
    facet_wrap(~researcher, nrow = 2) +
    labs(x = "Mail de rappel", y = "Pr(présent)")

Plusieurs manières de tester des hypothèses avec brms. La fonction hypothesis() calcule un evidence ratio (équivalent au Bayes factor). Lorsque l'hypothèse testée est une hypothèse ponctuelle (on teste une valeur précise du paramètre, e.g., \( \theta = 0 \)), cet evidence ratio est approximé via la méthode de Savage-Dickey. Cette méthode consiste simplement à comparer la densité du point testé accordée par le prior à la densité accordée par la distribution a posteriori.

(hyp1 <- hypothesis(mod1, "reminder = 0") )

Hypothesis Tests for class b:
      Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (reminder) = 0     1.56      0.16     1.22     1.88          0         0    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

1 / hyp1$hypothesis$Evid.Ratio

[1] 7041.326

plot(hyp1, plot = FALSE, theme = theme_bw(base_size = 20, base_family = "Open Sans") )[[1]] +
  geom_vline(xintercept = 0, linetype = 2) +
  coord_cartesian(xlim = c(-5, 5) )

data.frame(prior = hyp1$prior_samples$H1, posterior = hyp1$samples$H1) %>%
    gather(type, value) %>%
    mutate(type = factor(type, levels = c("prior", "posterior") ) ) %>%
    ggplot(aes(x = value) ) +
    geom_histogram(bins = 50, alpha = 0.8, col = "white", fill = "steelblue") +
    geom_vline(xintercept = 0, lty = 2, size = 1) +
    facet_wrap(~type, scales = "free") +
    labs(x = expression(beta[reminder]), y = "Nombre d'échantillons")

Une deuxième solution consiste à étendre l'approche par comparaison de modèles. Tester une hypothèse revient à comparer deux modèles : un modèle avec l'effet d'intérêt et un modèle sans l'effet d'intérêt.

prior2 <- c(
    prior(normal(0, 10), class = Intercept, coef = ""),
    prior(cauchy(0, 10), class = sd),
    prior(lkj(2), class = cor) )

mod2 <- brm(presence | trials(total) ~ 1 + reminder + (1 + reminder | researcher), 
    family = binomial(link = "logit"),
    prior = prior1,
    data = data,
    # this line is important for bridgesampling
    save_all_pars = TRUE,
    warmup = 2000, iter = 1e4,
    cores = parallel::detectCores(), backend = "cmdstanr",
    control = list(adapt_delta = 0.95) )

mod3 <- brm(presence | trials(total) ~ 1 + (1 + reminder | researcher), 
    family = binomial(link = "logit"),
    prior = prior2,
    data = data,
    save_all_pars = TRUE,
    warmup = 2000, iter = 1e4,
    cores = parallel::detectCores(), backend = "cmdstanr",
    control = list(adapt_delta = 0.95) )

On peut ensuite comparer la vraisemblance marginale de ces modèles, c'est à dire calculer un Bayes factor. Le package brms propose la méthode bayes_factor() qui repose sur une approximation de la vraissemblance marginale via le package bridgesampling (Gronau & Singmann, 2017).

bayes_factor(mod2, mod3)

Estimated Bayes factor in favor of mod2 over mod3: 2712061.97100

On peut également s'intéresser aux capacités de prédiction de ces deux modèles et les comparer en utilisant des critères d'information. La fonction waic() calcule le Widely Applicable Information Criterion (cf. Cours n°07).

waic(mod2, mod3, compare = FALSE)

Output of model 'mod2':

Computed from 32000 by 20 log-likelihood matrix

          Estimate  SE
elpd_waic    -59.6 2.2
p_waic        10.5 1.2
waic         119.2 4.4

15 (75.0%) p_waic estimates greater than 0.4. We recommend trying loo instead. 

Output of model 'mod3':

Computed from 32000 by 20 log-likelihood matrix

          Estimate  SE
elpd_waic    -58.9 1.6
p_waic        10.1 0.4
waic         117.9 3.2

19 (95.0%) p_waic estimates greater than 0.4. We recommend trying loo instead. 

Une autre manière d'examiner les capacités de prédiction d'un modèle est le posterior predictive checking (PPC). L'idée est simple : il s'agit de comparer les données observées à des données simulées à partir de la distribution a posteriori. Une fois qu'on a une distribution a posteriori sur \( \theta \), on peut simuler des données à partir de la posterior predictive distribution :

\[ p(\widetilde{y} | y) = \int p(\widetilde{y} | \theta) p(\theta | y) d \theta \]

Si le modèle est un bon modèle, il devrait pouvoir générer des données qui ressemblent aux données qu'on a observées (e.g., Gabry, Simpson, Vehtari, Betancourt, & Gelman, 2017).

On représente ci-dessous la distribution de nos données.

data %>%
    ggplot(aes(x = presence / total) ) +
    geom_density(fill = "grey20")

Cette procédure est implémentée dans brms via la méthode pp_check() qui permet de réaliser de nombreux checks. Par exemple, ci-dessous on compare les prédictions a posteriori (n = 100) aux données observées.

pp_check(mod2, nsamples = 1e2)

pp_check(mod2, nsamples = 1e3, type = "stat_2d")

En fittant des modèles un peu compliqués, il se peut que vous obteniez des messages d'avertissement du genre “There were x divergent transitions after warmup”. Dans cette situation, on peut ajuster le comportement de Stan directement dans un appel de la fonction brm() en utilisant l'argument control.

mod2 <- brm(
    presence | trials(total) ~ 1 + reminder + (1 + reminder|researcher), 
    family = binomial(link = "logit"),
    prior = prior2,
    data = data,
    warmup = 2000, iter = 1e4,
    cores = parallel::detectCores(), # using all availables cores
    control = list(adapt_delta = 0.95) # adjusting the delta step size
    )

On peut par exemple augmenter le pas de l'algorithme, via adapt_delta (par défaut fixé à 0.8), ce qui ralentira probablement l'échantillonnage mais améliorera la validité des échantillons obtenus. Plus généralement, soyez attentifs aux messages d'erreur et d'avertissement générés par brms.

Une liste de blogposts qui utilisent brms : https://paul-buerkner.github.io/blog/brms-blogposts/

L'article d'introduction du package brms (Bürkner, 2017) et la version “advanced” (Bürkner, 2018).

Un tutorial sur les modèles de régression logistique ordinale (Bürkner & Vuorre 2018).

Un article tutorial d'introduction aux modèles multi-niveaux avec brms (Nalborczyk et al., 2019, preprint) et un autre (Vasishth et al., 2018).

Une méta-analyse est simplement une analyse d'analyses. Il s'agit d'un modèle linéaire (presque) comme les autres, sauf que nos observations sont (généralement) des données déjà résumées par une taille d'effet (ou pas). On peut traiter ces tailles d'effets comme des observations, avec une variance connue.

On distingue deux classes de modèles.

• Méta-analyse à effets fixes (constant effects) : on considère que la taille d'effet estimée par toutes les études est la même
• Méta-analyse à effets aléatoires (varying effects) : on modélise la variabilité et les sources de dépendance dans les données

On ne considère que le deuxième type de modèle, en notant (cf. Cours n°08) qu'un modèle à effet fixe peut être considéré comme un modèle à effet aléatoire dont on a fixé \( \tau = 0 \).

Le jeu de données ci-dessous recense le résultat de 32 expériences visant à évaluer l'effet des contraintes bioméchaniques sur l'évaluation des distances (Molto, Nalborczyk, Palluel-Germain & Morgado, 2019).

d <- read.csv("data/meta.csv")
head(d, 15)

                              study experiment           yi         vi
1                   Costello (2015)          1  0.735209366 0.29846545
2                   Costello (2015)          2  0.925134052 0.04796738
3           Durgin & Russell (2008)          1 -0.124624069 0.17928177
4         Hutchison & Loomis (2006)          1 -0.128328394 0.10001576
5         Hutchison & Loomis (2006)          2  0.350487089 0.05963093
6             Kirsch & Kunde (2012)          1  0.609278390 0.31866670
7            Kirsch & Kunde (2013a)          1  0.511022255 0.16231426
8            Kirsch & Kunde (2013a)          2  0.407867633 0.04216585
9            Kirsch & Kunde (2013b)          1  0.229159882 0.32082216
10           Kirsch & Kunde (2013b)          2 -0.223592130 0.08876431
11           Kirsch & Kunde (2013b)          3  0.145865089 0.10038189
12 Lessard, Linkenauger & P. (2009)          1  0.172318408 0.09078810
13             Morgado & al. (2013)          1  0.002778791 0.13837508
14         Osiurak & Morgado (2012)          1  0.460865217 0.22233808
15  Proffitt, S, B & Epstein (2003)          1  0.351381281 0.10386001

On peut écrire un premier modèle de la manière suivante.

\[ \begin{aligned} y_{j} &\sim \mathrm{Normal}(\mu_{j}, \sigma_{j}) \\ \mu_{j} &= \alpha + \alpha_{study[j]} \\ \alpha_{study[j]} &\sim \mathrm{Normal}(0, \tau_{s}) \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \tau_{s} &\sim \mathrm{HalfCauchy}(0, 1) \\ \end{aligned} \]

Où \( \sigma_{j}^2 = v_{j} \) est la variance de l'effet dans la \( j \)ème étude et \( \alpha \) la taille d'effet dans la population. L'index \( \alpha_{study[j]} \) indique la taille d'effet moyenne dans l'étude \( j \). En plus de la variance dans l'échantillon \( \sigma_{j}^2 \) (qui est connue), on estime aussi la variabilité des effets entre les études \( \text{Var}(\alpha_{study}) = \tau_{s}^{2} \) (niveau 2).

Le modèle précédent peut être étendu à plus de deux niveaux pour prendre en compte les structures de dépendance dans le jeu de données. En effet, chaque étude (chaque papier publié) contenait plusieurs expériences. On pourrait s'attendre à ce que les expériences d'une même étude se ressemblent plus entre elles…

On peut écrire ce modèle sur trois niveaux, comme ci-dessous.

\[ \begin{aligned} y_{ij} &\sim \mathrm{Normal}(\mu_{ij}, \sigma_{ij}) \\ \mu_{ij} &= \alpha + \alpha_{study[j]} + \alpha_{experiment[ij]} \\ \alpha_{study[j]} &\sim \mathrm{Normal}(0, \tau_{s}) \\ \alpha_{experiment[ij]} &\sim \mathrm{Normal}(0, \tau_{e}) \\ \alpha &\sim \mathrm{Normal}(0, 1) \\ \tau_{e}, \tau_{s} &\sim \mathrm{HalfCauchy}(0, 1) \\ \end{aligned} \]

On estime maintenant, en plus de la variabilité \( \sigma_{ij} \), deux autres sources de variation : la variation des effets entre différentes expériences issues d'une même étude \( \text{Var}(\alpha_{experiment}) = \tau_{e}^{2} \) (niveau 2) et la variation entre différentes études \( \text{Var}(\alpha_{study}) = \tau_{s}^{2} \) (niveau 3).

Ce modèle se fit facilement avec brms.

prior4 <- c(
    prior(normal(0, 1), coef = intercept),
    prior(cauchy(0, 1), class = sd)
    )

mod4 <- brm(
    yi | se(sqrt(vi) ) ~ 0 + intercept + (1 | study) + (1 | experiment),
    data = d,
    prior = prior4,
    save_all_pars = TRUE,
    warmup = 2000, iter = 1e4,
    cores = parallel::detectCores(),
    control = list(adapt_delta = 0.99),
    backend = "cmdstanr"
    )

summary(mod4)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: yi | se(sqrt(vi)) ~ 0 + intercept + (1 | study) + (1 | experiment) 
   Data: d (Number of observations: 32) 
  Draws: 4 chains, each with iter = 8000; warmup = 0; thin = 1;
         total post-warmup draws = 32000

Group-Level Effects: 
~experiment (Number of levels: 4) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.36      0.26     0.06     1.03 1.00    10322    11288

~study (Number of levels: 19) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.22      0.10     0.04     0.43 1.00    10081    10609

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
intercept     0.42      0.23    -0.04     0.89 1.00    13160    13273

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.00      0.00     0.00     0.00   NA       NA       NA

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

mod4 %>%
  plot(
    pars = c("^b_", "^sd_"),
    combo = c("dens_overlay", "trace"),
    theme = theme_bw(base_size = 16, base_family = "Open Sans")
    )

La statistique bayésienne est une approche générale de l'estimation de paramètres. Cette approche utilise la théorie des probabilités pour quantifier l'incertitude vis à vis de la valeur des paramètres de modèles statistiques.

Ces modèles sont composés de différents blocs (e.g., fonction de vraisemblance, priors, modèle linéaire ou non-linéaire) qui sont modifiables à souhait. Ce qu'on appelle classiquement “conditions d'application” sont simplement les conséquences des choix de modélisation réalisés par l'utilisateur. Autrement dit, c'est l'utilisateur qui choisit (et ne subit pas) les conditions d'application.

Nous avons vu que le modèle de régression linéaire est un modèle très flexible qui permet de décrire, via la modification de la fonction de vraisemblance et via l'introduction de fonctions de lien, des relations complexes (e.g., non-linéaires) entre variable prédite et variables prédictrices. Ces modèles peuvent gagner en précision par la prise en compte de la variabilité et des structures présentes dans les données (cf. modèles multi-niveaux).

Le package brms est un véritable couteau suisse de l'analyse statistique bayésienne en R. Il permet de fitter presque n'importe quel type de modèle de régression. Cela comprend tous les modèles que nous avons vu en cours, mais également bien d'autres. Entre autres, des modèles multivariés (i.e., avec plusieurs outcomes), des modèles “distributionnels” (e.g., pour prédire des différence de variance), des generalised additive models, des procesus Gaussiens, des modèles issus de la théorie de détection du signal, des mixture models, des modèles de diffusion, des modèles non-linéaires…

N'hésitez pas à me contacter pour plus d'informations sur ces modèles ou si vous avez des questions par rapport à vos propres données. Vous pouvez aussi contacter le créateur du package brms, très actif en ligne (voir son site). Voir aussi le forum Stan.

Ce jeu de données recense des données concernant 2000 élèves dans 100 écoles différentes. L'outcome principal est la popularité de l'élève, évaluée sur une échelle de 1 à 10, et estimée en utilisant une procédure sociométrique (i.e., on a demandé aux élèves de se noter mutuellement). Ces élèves étaient également notés par leurs professeurs (colonne teachpop), sur une échelle de 1 à 7. On dispose comme prédicteurs du genre de l'élève (boy = 0, girl = 1) et de l'expérience du professeur (texp, en années).

d <- read.csv("data/popular.csv")
head(d, 10)

   pupil school popular  sex texp teachpop
1      1      1       8 girl   24        7
2      2      1       7  boy   24        7
3      3      1       7 girl   24        6
4      4      1       9 girl   24        6
5      5      1       8 girl   24        7
6      6      1       7  boy   24        7
7      7      1       7  boy   24        7
8      8      1       7  boy   24        7
9      9      1       7  boy   24        7
10    10      1       8  boy   24        6

À vous d'explorer ce jeu de données, de fitter quelques modèles (avec brms) pour essayer de comprendre quels sont les facteurs qui expliquent (permettent de prédire) la popularité d'un élève…

d %>%
    ggplot(aes(x = popular) ) +
    geom_histogram() +
    facet_wrap(~sex) +
    scale_x_continuous(breaks = 1:10, limits = c(1, 10) )

d %>%
    ggplot(aes(x = texp, y = popular) ) +
    geom_point(alpha = 0.2) +
    geom_smooth(method = "lm", colour = "black") +
    facet_wrap(~sex)

library(magrittr)

d %<>%
    mutate(
        # contrast-coding gender
        sex = ifelse(sex == "boy", -0.5, 0.5),
        # centering and standardising teacher experience
        texp = scale(texp) %>% as.numeric
        )

prior5 <- c(
    prior(normal(5, 2.5), class = Intercept),
    prior(cauchy(0, 10), class = sd),
    prior(cauchy(0, 10), class = sigma)
    )

mod5 <- brm(
    popular ~ 1 + (1 | school),
    data = d,
    prior = prior5,
    save_all_pars = TRUE,
    warmup = 2000, iter = 1e4,
    cores = parallel::detectCores()
    )

prior6 <- c(
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 1), class = b),
    prior(cauchy(0, 1), class = sd),
    prior(cauchy(0, 10), class = sigma)
    )

mod6 <- brm(
    popular ~ 1 + texp + (1 | school),
    data = d,
    prior = prior6,
    save_all_pars = TRUE,
    warmup = 2000, iter = 1e4,
    cores = parallel::detectCores()
    )

prior7 <- c(
    prior(normal(0, 1), class = Intercept),
    prior(normal(0, 1), class = b),
    prior(cauchy(0, 1), class = sd),
    prior(cauchy(0, 10), class = sigma),
    prior(lkj(2), class = cor)
    )

mod7 <- brm(
    popular ~ 1 + sex + texp + (1 + sex | school),
    data = d,
    prior = prior7,
    save_all_pars = TRUE,
    warmup = 2000, iter = 1e4,
    cores = parallel::detectCores()
    )

mod8 <- brm(
    popular ~ 1 + sex + texp + sex:texp + (1 + sex | school),
    data = d,
    prior = prior7,
    save_all_pars = TRUE,
    warmup = 2000, iter = 1e4,
    cores = parallel::detectCores()
    )

# calcul du WAIC et ajout du WAIC à chaque modèle
mod5 <- add_criterion(mod5, "waic")
mod6 <- add_criterion(mod6, "waic")
mod7 <- add_criterion(mod7, "waic")
mod8 <- add_criterion(mod8, "waic")

# comparaison des WAIC de chaque modèle
w <- loo_compare(mod5, mod6, mod7, mod8, criterion = "waic")
print(w, simplify = FALSE)

     elpd_diff se_diff elpd_waic se_elpd_waic p_waic  se_p_waic waic    se_waic
mod8     0.0       0.0 -1990.3      32.5        161.1     5.2    3980.6    65.1
mod7    -1.8       2.0 -1992.1      32.7        163.9     5.3    3984.2    65.4
mod6  -448.7      25.8 -2439.0      30.2         93.2     2.8    4877.9    60.5
mod5  -449.6      25.9 -2439.9      30.3         95.1     2.8    4879.8    60.5

Les prédictions du modèle ne coïncident pas exactement avec les données car celles-ci sont catégorielles. Les élèves étaient notés sur une échelle discrète allant de 1 à 10 (i.e., un élève ne pouvait pas avoir une note de 3.456). Ce type de données peut-être approximée par une distribution normale (comme nous l'avons fait) mais ce choix n'est pas optimal en termes de prédiction…

pp_check(mod8, nsamples = 1e2)

On pourrait choisir un modèle qui se rapproche du processus de génération des données. C'est le cas du modèle de régression logistique ordinale, ou ordered categorical model. Ce modèle est une sorte de généralisation à plus de 2 catégories du modèle de régression logistique vu au Cours n°06 (voir ce blogpost pour plus de détails, ou McElreath, 2015, chapitre 11), sauf que les catégories sont ordonnées.

\[ \begin{aligned} \text{pop}_{i} &\sim \mathrm{Categorical}(\mathbf{p}) \\ \text{logit}(p_{k}) &= \alpha_{k} \\ \alpha_{k} &\sim \mathrm{Normal}(0, 10) \end{aligned} \]

Où la distribution \( \mathrm{Categorical} \) est une distribution discrète qui prend un vecteur de probabilités \( \mathbf{p} = \{p_{1}, p_{2}, p_{3}, p_{4}, p_{5}, p_{6}, p_{7}, p_{8}, p_{9}\} \) qui correspondent aux probabilités cumulées de chaque réponse (entre 1 et 10, 10 ayant une probabilité cumulée de 1).

On définit une série de \( N - 1 \) intercepts \( \mathbf{p} = \{p_{1}, p_{2}, p_{3}, p_{4}, p_{5}, p_{6}, p_{7}, p_{8}, p_{9}\} \) sur le logarithme de la cote cumulée (log-cumulative-odds).

\[ \text{logit}(p_{k}) = \log \frac{\Pr(y_{i} \leq k)}{1 - \Pr(y_{i} \leq k)} = \alpha_{k} \]

La vraisemblance de l'observation \( k \) (e.g., pop = 3) est donnée par soustraction des proportions cumulées. Cette vraisemblance est représentée par les barres verticales sur le graphique ci-dessous.

\[ p_{k} = \Pr(y_{i} = k) = \Pr(y_{i} \leq k) - \Pr(y_{i} \leq k -1) \]

NB : Ce modèle peut prendre plusieurs (une dizaine) heures…

mod9 <- brm(
    popular ~ 1 + sex + texp + sex:texp + (1 | school),
    data = d,
    prior = prior6,
    warmup = 2000, iter = 5000,
    chains = 4, cores = parallel::detectCores(),
    file = "models/mod9", backend = "cmdstanr"
    )

prior10 <- c(
    brms::prior(normal(0, 10), class = Intercept),
    brms::prior(normal(0, 10), class = b),
    brms::prior(cauchy(0, 10), class = sd)
    )

mod10 <- brm(
    popular ~ 1 + sex + texp + sex:texp + (1 | school),
    data = d,
    family = cumulative(link = "logit"),
    prior = prior10,
    chains = 4, cores = parallel::detectCores(),
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "models/mod10", backend = "cmdstanr"
    )

waic(mod9, mod10, compare = FALSE)

Output of model 'mod9':

Computed from 12000 by 2000 log-likelihood matrix

          Estimate   SE
elpd_waic  -2086.5 31.3
p_waic        96.7  3.0
waic        4173.0 62.6

7 (0.4%) p_waic estimates greater than 0.4. We recommend trying loo instead. 

Output of model 'mod10':

Computed from 4000 by 2000 log-likelihood matrix

          Estimate   SE
elpd_waic  -2074.7 32.4
p_waic       106.5  2.6
waic        4149.4 64.9

3 (0.1%) p_waic estimates greater than 0.4. We recommend trying loo instead. 

pp_check(mod10, nsamples = 1e2, type = "bars", prob = 0.95, freq = FALSE) +
  scale_x_continuous(breaks = 1:9) +
  labs(x = "Popularité", y = "Proportion")